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ABSTRACT 

In this work we investigate the multivariate statistical description of the matter distribution 
in the nonlinear regime. We introduce the multivariate Edgeworth expansion of the lognor- 
mal distribution to model the cosmological matter field. Such a technique could be useful 
to generate and reconstruct three-dimensional nonlinear cosmological density fields with the 
information of higher order correlation functions. We explicitly calculate the expansion up 
to third order in perturbation theory making use of the multivariate Hermite polynomials up 
to sixth order. The probability distribution function for the matter field includes at this level 
the two-point, the three-point and the four-point correlation functions. We use the hierarchi- 
cal model to formulate the higher order correlation functions based on combinations of the 
two-point correlation function. This permits us to find compact expressions for the skewness 
and kurtosis terms of the expanded lognormal field which can be efficiently computed. The 
method is, however, flexible to incorporate arbitrary higher order correlation functions which 
have analytical expressions. The applications of such a technique can be especially useful to 
perform weak-lensing or neutral hydrogen 21 cm line tomography, as well as to directly use 
the galaxy distribution or the Lyman-alpha forest to study structure formation. 
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1 INTRODUCTION 

The cosmological matter distribution encodes the information of 
the origin of our Universe and the processes which lead to structure 
formation. A precise understanding and modeling of its statistics is 
crucial to extract the cosmological information from observational 
data and to ultimately understand cosmic evolution. 

The Universe is being scrutinized with unprecedented accu- 
racy. Many excellent astronomical surveys have been launched in 
the recent past and ongoing and upcoming projects are on the way 
to perform the most ambitious map of the Universe up-to-date. 

Some of the most successful lo w-redshift galaxy ca talogs are 
the 2dF Galaxy Redshift SurveyR dColless et all l2003h a nd the 
Sloan Digital Sky Survey (SDSS0 1 Abazaiian et al.ll2009h - Deei 



scope Legacy Survey (CFHTLSfl dHoekstra et all l2006h . the 
VISTA (Visible and Infrared Survey Telescope for Astron- 
omyQ and the PANoramic Sur vey Telescope And Rapid 
Respo nse System (Pan-STARRSfl dKaiser & Pan-STARRS Tea 
20020. and the planned Dark E n ergy Survey (DES 
I The Dark Energy Survey C ollaboration 2005) or the Joint 



surveys like the Bary on Oscillation Spectros copic Survey (BOSS 
dSchlegel et aT]|2009h . the DEEP2 Survey Fl dDavis et al.ll2003h and 



the VIMOS VLT Deep Survey (VVDS)Q(Le Fe vre et al.l20oi) are 
being run. 

On the other hand the Canada-France-Hawaii Tele- 



Dark Energy Mission (JDEM13 from NA SA-DOE and Dark 
UNiverse Explorer (DUNE) from CNES dCrotts et alj 120051; 
iRefregier et al J 120061) will provide the first weak lensing surveys 
covering very large sky areas and depth. 

Also the Lyman alpha forest will be a useful observable to 
study cosmology using for instance the BOSS survey (the matter 
power-spectrum has a lready been measured with the SDSS, see 
iMcDonald et alj2005h . 

The neutral hydrogen 21 -cm line will provide a new astronom- 
ical window to study structure formation in the Universe. Some of 
the most notable projects are the Giant Metre- wave Radio Tele- 
scope (GMRtEI dPen et aljUooj) . the Precision Array to Probe 
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Epoch of Reionization (PAPERf 12 ] ( P arsons et ai1l2010h . the LOw 
Frequency ARray (LOFARpl ( Falcke et alj2007h an d the Murchi- 
son Widefield Array (MWAtj [ Lonsdale et alj|2009h . 

In summary, an avalanche of astronomical data is being col- 
lected to study its structure and history based on different observ- 
ables. In order to extract valuable cosmological information not 
only a careful modeling of the systematics of the observation pro- 
cess and the nature of the observable is required, but also a precise 
modeling of the underlying signal. We focus here on the cosmolog- 
ical matter density field. 

Different approaches can be found in the literature to recon- 
struct the large-scale structure. Geometrical reconstruction meth- 
ods try to approximately capture the higher order statistics beyond 
the two-point correlation function in an effective way through a 
geometry-based prescription to form structures from a point source 
distribution. The salient and pervasive foamlike pattern of the cos- 
mic web has led to develop method s like the Voronoi or De- 
launay tessellations (see f o r example van de Weygaert et alj 20091 ; 
lAraeon-Calvo et alj [20071 ; ISchaap & van de Weygaertl l200ol) . On 
the other hand one can find reconstruction methods based on a sta- 
tistical approach. The advantage of the statistical methods with re- 
spect to the geometrical ones is that one can clearly specify the 
assumptions made on the matter field and the observable in form 
of probability distribution functions. The statistical methods could 
be more suitable to extract statistical quantities like the power- 
spectrum. This has been succesfully done in the Cosmic Microwave 
Background (CMB) for which the fluctuatio ns can be assumed to 
be Gaussian distributed dEriksen et alj|2007h . The disadvantage is 
that the level of complexity that such methods can achieve has 
always been limited as the formulation of the probability distri- 
bution functions and its applications to reconstruction methods is 
relatively complex and computer-intensive. Indeed, any complex 
realistic formulation seemed to be untreatable. For a long time 
the Wiener-filter, also-called least-squares filter, has been the only 
available method in Astronomy to incorporate the stati stical infor- 
mation in the reconstruction method (see for example Bunn et all 
1994 IZaroubi et al.ll 1991 iFisher et al.ll 19951; IWebster et al.ll 19971 : 
Zaroubi et a l. 1999; Erdogdu et al. 2004, 2006). Less attention was 
paid in the astrophysics commun ity to the nonlinear version of the 
least-squares filter proposed by Tar antola & Valettel dl982l) . Here 
the data model which relates the measurements to the seeked sig- 
nal is extended to be nonlinear. Probably the first group apply- 
ing this method in an astrophysical context was Chen et al. ( 1998) 
to map the interste llar absorption structures in the galactic plane. 
Picho n et alj j200 lh proposed to use this nonlinear reconstruction 
scheme to recover the cosmic density and velocity field traced 
by the Lyman alpha forest. The drawback of iTarantola & Valettel 
(1982)'s approach is that it requires both a Gaussian prior (with a 
nonlinear transformation) and a Gaussian likelihood for the dis- 
tribution of the observable. This is a too crude assumption for 
many observables, like for instance a galaxy distribution. Recently 
the Poisson-lognormal (and Gaussian-lognormal) model was pro- 
posed in a Bayesian framework to recover the cosmic density field 
dKitaura et aflWol) . In this study it was shown that the lognor- 
mal prior is in good agreement with the underlying matter field 
extracted from N-body simulations in the large overdense regions 
(> 10 3 ), but fails to fit the matter statistics in the underdense 



12 http://astro.berkeley.edu/ dbacker/eor/ 
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14 http://www.MWAtelescope.org/ 



regions. One of the advantages of this method with respect to 
the nonlinear least-squares approach is that it can deal with non- 
Gaussian likelihoods. Another important point of the Bayesian ap- 
proach proposed in dKitaura etatll201 Oh is that it can be easily 
extended to sample full posterior distributions (see the wo rks by 
IJasche & Kitaur j|2010l ; [jasche et al.ll20ld : lKitaura et alfcoioh . 

Nevertheless, non of the above mentioned statistical methods 
includes any information beyond the two-point correlation func- 
tion. As gravitational clustering forms nonlinear structures the Uni- 
verse becomes inhomogenous and complex patterns arise which en- 
code high order statistics. 

The purpose of this work is to extend the multivariate char- 
acterization of matter beyond the two-point correlation function to 
incorporate higher order statistics. We therefore relax the lognor- 
mal assumption and introduce the multivariate Edgeworth expan- 
sion which leads to additional terms describing the skewness and 
kurtosis of the field with respect to the lognormal distribution func- 
tion. This wor k is based on the univar i ate Ed g eworth expansion 
introd u ced by Scherrer & Bertsch inger (1991); Juszk iewicz et al.1 
dl995h : iBernardeau & Kofmanl dl995h and IColombil dl994h . The 
Edgeworth expansion we find deviates from the trivial general- 
ization of the univariate c ase to the multivariate case. We use 
the hierarchical model (see iFrv&Peeblesll 19781 ; lFrvlll984 1 19861 : 
Balian & Schaef fei | |l989h to formulate the three-point and four- 
point correlation functions which permits us to find particular ex- 
pressions for the skewness and kurtosis terms. The expressions 
we find are compact due to the symmetries introduced by the 
hierarchical model and can be computed by means of convolu- 
tions with fast Fourier transforms (fft's). As the hierarchical model 
may fail at certain scale s and regimes ( Suto & Matsubara 1 19941 : 
iMatsubara & Sutoll 19941) this work could be extended incorporat- 
ing more complex higher order correlation functions which in- 
clude galaxy biasing or redshift distortions (see the works by 
I Scoccimarro et alj 1 19981; iTaruva & Sodal 1 19991 ; IMatsubara! 120031 : 
IZheng|2004lMatsubaral2008l) and t o perform topological and mor- 
phological studies (see for examp le Matsub ara & Yokovama|[l99r3 : 
iGott et al.ll2008l ; |james e"tai]|2009h . 

We believe that the method introduced in this paper can be 
very useful to study cosmological structures in the range between 
the quasi-nonlinear and the nonlinear regime. It could be interest- 
ing to apply higher order statistics to galaxy redshift surveys, to 
weak-lensing surveys, to the Lyman alpha forest or to the 21 cm 
line. We would like to warn the reader that this work is still in a 
development phase as higher order correlation models need to be 
tested and many numerical studies still have to be done. This is the 
first of a series of works in which we will analyze the statistical 
description of gravitational clustering. 

This paper is structured as follows. In the next section the 
gravitational clustering statistics will be analyzed in great detail 
(section [2}. We will start reviewing the work done so far for the 
univariate matter distribution (section 12.2b and then present the 
multivariate case (section [23). This will lead us to a multivariate 
Edgeworth expansion of the lognormal field up to third order in 
perturbation theory including two-point, three-point and four-point 
correlation functions. Then we will present the hierarchical model 
(section |2"l4T > and use the expression for the three-point correlation 
function to calculate the skewness and kurtosis terms in the Edge- 
worth expansion (section |2~5l l. A detailed calculation can be found 
in the appendix. Finally we will present the summary and conclu- 
sions of this work. 



© 0000 RAS, MNRAS 000, 000-000 



Non-Gaussian field statistics 3 



2 GRAVITATIONAL CLUSTERING FIELD STATISTICS 

In this section we will study the matter field statistics produced by 
gravitational clustering. We start with a physical motivation fol- 
lowed by the review of the univariate non-Gaussian statistics. Then 
we introduce the multivariate Edgeworth expansion of the Lognor- 
mal field. Finally we present the hierarchical model and calculate 
the skewness and kurtosis terms of the Edgeworth expansion. 



2.1 Physical motivation 

Let us divide the Universe into N c cells and assign to each cell i 
a position r,, a matter density pi and a peculiar velocity Vi. The 
continuity equation relates the evolution of the matter content in 
the Universe to its peculiar velocity field: 







dp 1 , ,„ 



(1) 



with t being the cosmic time, a the scale factor, r the set of posi- 
tions ({ri , . . . , rjv c }), p the matter density field ({pi , . . . , pjv c }), 
v the peculiar velocity field ({vi , . . . , vn c }) and d/dt = d/dt + 
1/a (v ■ Vr)p(r) the total derivative. 

We can follow matter particles until they start crossing-over 
(in this regime particles can have different peculiar velocities at the 
same position) and form caustics. Before this occurs we can write 
the formal solution to Eq. (fTJ as: 



(p)e s , s 



dt-Vr ■ v, 

a 



with s being the logarithm of the normalized density: 
s = lnp- ln(p) = m(p/(p}) = ln(l + 



(2) 



(3) 



and the matter overdensity field given by: Sm = p/{p) — L The 
ensemble averages are used at this stage only to denote the mean of 
the variable. Note that the field s does not have zero mean, but is 
given by: p s = (s) — (In p) — ln(p). It is convenient to define a 
field <I> with zero mean ((<&) = 0): 



$ = In p - (In p) = s - n s 



(4) 



Assuming that $ is a Gaussian ra ndom field leads to a log- 
normal distributed density field (see IColes & Joneslll99lh . Note 
however, that Lagrangian perturbation theory -which is known 
to give a good approximation of gravitational cluster i ng until 
shell crossing starts (see e.g. iBuchert & Ehlerj 1 19931 ; [Buchert 
1 19941 : iBouchet et alj Q995)- dev iates from the l ognormal distri- 
buti on already in the linear IZel'dovichl d 1 97dh approximation 
(see Padmanabhan & Subramanianl 1 1993c [Bernardea u & Kofmanl 
1 1995 ). Furthermore, after structures start to virialize the peculiar 



velocity field will be strongly modified a nd Lagrangian p ertur- 
bation theory will start to fail dramatica lly. IColombil d 19941) sug- 
ges ted to use the formalism develo ped bv ljuszkiewicz et alj i 19951) 
and lBernardeau & Kofm an ( 1995) to study the departures from the 
lognormal distribution function including higher order correlation 
functions in the univariate matter distri bution (see also the w ork on 
a generalized lognormal distribution bv lSzapudi et alfe OQO), 



2.2 Univariate case 

In this subsection we will revise the matter statistics for the 
one-dimensional probability distribution function as developed 



in the work s by Juszkiewi cz et al. | d 19951) : iBernardeau & Kofmanl 
dl995h and lColombil d 19941) . F or a general overview on a symp- 
totic statistical techniques see Berkowitz & Garneil dl970l) and 
iBarndorff-Nielsen & Coxl Il989). Other univariate matter field dis- 
tribution functions have been proposed. They are however not triv- 
ially extendable to the multivariate case. Either they do not in- 
clude higher order correlations, but are extracted from fitting the 
univariate matter distribution based on numerical N-body simula- 
tions (e.g. Miralda-Es cude et alj|2000h , or a distribution function 
function is expanded u sing the variance as an univariate parameter 
dGaztanaga et alj |2000). For this reason we will consider only the 
above mentioned approach based on the expansion of the lognor- 
mal distribution function. 

Let us define the quantity v with zero mean and unity variance: 



(5) 



with a 2 — ("I> 2 ) being the variance of $ (the relation between the 
variance of $ and the variance of the matter overdensity 8m is de- 
rived in appendix |Bj. Higher order moments of v can be found by 
calculating the ensemble average of powers of v over the probabil- 
ity distribution function P(u): 



Pn 



dvP{y) v n = (O, 



(6) 



with n being the order of the moment. Please note that the moment 
/i refers to the variable v and not to the variable s. The moment 
generating function is given by: 



Mu{t) = Y.^- = / ^p^y tv = (e tv ) 



(7) 



Subsequent derivatives of M v (t) at the origin t = yield the mo- 
ments: 



d n M„{t) 



At" 



The cumulant generating function is given by: 



C (*) = I>»-p 



with n n being the cumulants or connected moments: 

K n = {is n ) c . 



(8) 



(9) 



(10) 



The cumulants can be obtained from the relation between the 
moment and cumulan t generating functions (see for example 
IBernardeau et alj|2002l) : 



M„(t) = exp(C(t)), 



(11) 



or equivalently: C(t) = In (A4 v (t)). Hence, the cumulants can be 
calculated by: 



d n In (M v (t)) 



dt" 



(12) 



It is however, more convenient to use expression i ll It to relate the 
moments to the cumulants: 



n=0 \n=l / 



(13) 
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This equation yields for the first order moments 
1 

Kl = 

2 , 
H2 + Ki = 1 

K3 + 3K2K1 + «! = 

K4 + 4K3K1 



/'•() 
Ml 
M-> 
A»3 
/i4 

Ms 



3^2 + 6fC2 Kl + K l — K 4 + 3 



K5 + 5K4K1 + IOK3K2 + IOK3K1 + 15K2K1 + lOft^^l 



Pd = K6 + 6K5K1 + 15ft4K2 + 15^4^1 + 10k 3 + 6OK3K2K1 
+ 2OK3K1 + 15fv2 + 45^2^1 + 15k2Ki + Ki 
= K6 + 15«4 + IOK3 + 15 . 

The probability distribution function of v can be obtained by 
inverting Eq. Q using the inve rse Laplace transform (see the re- 
view by Bernard eau et alj2002h : 



P{v) 



At 



2liy 



:exp{tv + C{t)) . 



(15) 



Thus, the generating function fully defines the probability distribu- 
tion function. When the departures from the Gaussian distribution 
function are small one can expand P(y) in a Gram-Charlier series: 



P(u) = G[y) 



with G(v) - — 

y ' v27T 

The Hermite polynomial of degree n can be calculated by 



1 2 H n 1 2 



This leads to the following first polynomials: 



hi(v 

h z (v 
hi(y 
hs(v 
h 6 {v 



1 

6v 2 + 3 
10p 3 + Ibv 



15. 



if I 7^ m; 
otherwise, 



; - 15i> 4 + 45z/ 2 
Using the orthogonality relation: 

&v G{v)hi(v)h m {v) = 

we can calculate the Gram-Charlier coefficients c; : 

C! = (-1) ! j &vP{v)hi{v) = (-lf(hi(u)) 

The latter expression yields: 

1 




— = — ^3 

/X4 — 3 = ^4 

—fX5 + 10^3 = — KB 



l + ^2pCi(-l)%(v) . (16) 

!=1 

= e~ "2" and hi (y) being the Hermite polynomials. 

(17) 



(18) 



(19) 



(20) 



(21) 



The Gram-Charlie r series may have poor convergence properties 
(seelCramej|l946[), and in so me cases even viol ently diverge (see 

( ' iBlinnikov & Moessned Il998h . For this reason ljuszkiewicz et al.l 
1 19951) suggested to model the departures from the Gaussian dis- 
tribution with an Edgeworth expansion which is a true asymp- 
totic expansion (for a full ex plicit expansion for arbitrary order see 
IBlinnikov & Moessnerl 1 199a . and references therein). The Edge- 
worth expansion consists of rearranging the terms in the Gram- 

■f Ki Charlier series based on collecting all the terms with the same clus- 
tering strength. To see how to do this let us recall how the different 
terms scale with a. 

In perturbation theory the matter field is expanded as a sum 
of terms with increasing order: $ = $1 + $2 + $3 + • • • , with 
$1 = 0(a 1 ) being the linear term, $2 = C(cr 2 ) being the 
quadratic term, $3 = 0(a 3 ) being the cubic term, etc.. The lin- 
ear term $1 is assumed to be Gaussian distributed, or equivalently 
in our case, the matter density field is assumed to be lognormal 
distributed at first order. To calculate the subsequent perturbation 
terms one nee ds to u se higher order correlation functions. As it was 
shown bv lFrvl dl984l) the cumulants of order higher than two scale 
as: Kyi = 0{a 2 "- 2 ). For this reason the following normalized cu- 
mulant has usually been introduced to calculate the Edgeworth ex- 
pansion: 



(22) 



with S3 = and S4 = ^ ' > ~&'' T ■ Using the latter expression 
one can group the terms with the same scaling power of a. The 
Edgeworth expansion until third order perturbation is then given 
by: 

P (4>) = % . P{v) = a-' J= exp (- {a ~y 



1 + [y Sihsia- 1 ®) 
+a 2 f iSi/u(<7 -1 *) + ^Ssheia- 1 ^)) + ... . (23) 



2.3 Multivariate case 

In this subsection we generalize the relations of the univariate mat- 
ter distribution to the multivariate case. Now <&, p and s are scalar 
fields: 



and 



$i = In f>i - (In p) 



Ms , 



(24) 



(25) 

= (ln(l + 6ui] ln(l + 6uj)) - (ln(l + 5 M i))(ln(l + S Mj )) , 

with S being the variance of the field ln(l + <5M) — Ps (see appendix 
[C]for the relation between the variance of 3? and the variance of the 
matter overdensity field Sm)- We introduce the field u which has 
zero mean and unity variance by definition: 



The n-dimensional moments are given by: 



pa — 15/^4 + 30 = Ke + Wk 3 . 



&i>P(v) v tl 



■Via)- 



(26) 



(27) 
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□ 



l^ijk 



f-^ijkl 



H'ij klm. 



f-^ijklmn 




^ i j k I. 





' + / + \ + 

• • • • • • • 



f^r'l ^ J ^k 



15 



10 



K'ijklm ^ijkl^m ^ijk^lm ^ijk^l^m ^ijf^kl^m ^ij^k^l^m- f^i^j f^k^l^m 



° 



15 




60 



20 




^ijklmn ^ijklm^n ^ijkl^mn f^ijkl^m^n ^ijk^lmn ^ijk^lm^n ^ijk^l^m^n f^ij^kl^mn 



45 



15 



Figure 1. Moments up to six point correlation statistics. The rectangular box stands for the ensemble average over the points, i.e. the moments. The points 
connected through lines represent the connected moments. In case of more than one equivalent element the number of elements is indicated. The circles mark 
the terms which are not automatically zero as they do not contain a single unconnected point (we deal here with a centered variable with zero mean). The 
three equivalent terms for the third order moment have been specified. Let us define a cumulant object as a set of connected points in a term of a moment. Two 
cumulant objects will be equivalent if they have the same number of points. The number of equivalent elements N c in a term is calculated by the factorial of 
the order of the moment n divided by the factorial of the number of single points 7V p o , the factorial of the number of equivalent cumulant objects iV bj in a 
term and the factorial of the number of points for each cumulant object T7 fc N p k : N c = nl /(N p o ! iV bj ! Ylk ^pk '•)• As an example the eighth term of the 
sixth order moment: KijK k iK mn has 15 equivalent elements: 15 = 6!/(3!2!2!2!). 



The multivariate moment generating function yields: 



Mv{t\ ■ ..t„) 



Sl---9n=0 



L \ • • • L n 



(28) 



Analogously to the univariate case subsequent derivatives of 
Mu(t) at the origin t — lead to the moments: 



d n Mu{ti...t n ) 



dti . . . dt n 



t 1 ...t n =o 



The cumulant generating function is given by: 



q 1 ...q n =0 



fH fin. 



(29) 



(30) 



with Ki 1 ...i n being the cumulants or connected moments: 

Ki!...i„ = {Vi-L ■ ■ ■ Vi n )c ■ (31) 

The moments are related to the cumulant generating functions by: 
Mv(ti ...*„) = exp(C(ii . . . t n )) . (32) 
Plugging in Eqs. (28) and d30t in the latter expression we obtain: 

<}l...q„=0 

exp ( E o (<--t) c fr-5)- (33) 

\91...8n=0 * / 

A way to spare the tedious calculations consists on looking at 
all combi nations of connections b etween points in a diagram (see 
Fig. [U and iBernardeau et al J 20021) . 
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Since the mean of Vi is zero {p,i = 0) we got rid off all the 
terms containing single connected points (/tj). The result for the 
first moments is listed below: 



Here are the results for the first couple of polynomials: 



/■'(> 


= 1 


Mi 


= Ki = 


/''J 




l^ijk 


— K ij k 







(34) 



E 

ji... J4 e[i,...,4] 



f-l'ijklm, — fcijklr, 



+ 



- T 

3!2 ^ 



f^ijklmn 



l 

4!2 



ji... J5 e[i,...,5] 

E 



3!3!2 



E 



3 1 ...j 6 6[l,...,6] 



+ 



2 3 



/ J J1---J6 *j 1 l j 2 l j3 l 34 '• 



K 



ii...j 6 e[i,...,6] 



where we have used the following identities: i = ii,j = i-2, k = 
i3,l = 44, m = is, n = ie , the Kroenecker delta: S K and the 
modified Levi-Civita tensor we introduce here: 



( 1) tei j 1 --- i j n 



(35) 



with N t being the number of transpositions. This tensor has the 
property of being always positive (including zero) since it is multi- 
plied by the factor (— 1) * which is positive for an even number of 
transpositions and negative for an odd number of transpositions, 
thus compensating for the negative sign coming from the Levi- 
Civita tensor. The number of equivalent objects is indicated to the 
lower right of the terms. To see how this number is calculated see 
the caption in Fig. (TJ. 

Let us study here the multivaria te Gram-Charlier series expan- 
sion (see lBerkowitz & Garnerlll970l) : 



P{u) = G{u) 



°° 1 

1+E E n^l-^lh 1 ) 1 ^!-^^) 



(36) 



with G{v) being a multivariate Gaussian distribution G(u) — 

1 _vH>_ 
-k=e 2 ■ 

V2tt 

The corresponding multivariate Hermite polynomials are cal- 
culated by (see lBerkowitz & Garnej[T97(]|) : 



dfi 1 . . . dv, 

from which the following recursive formula can be built: 



e-i» f », (37) 



» = (-l) w e = 



I i\n-l + 

(-1) e 2 hi 



(38) 



which we have used in our calculations. 



h {u) = 1 
hi (y) = Vi 
hij(y) = v l v j - Sfj 



(39) 



ViVjV k 



3132336[1,2,3] 
hijkl{v) = ViVjV k Vi 



- y 

)2 Z_^i 



2- 



rK 

£,*, Vi . Vi- Oj . ,■ . 



il...j 4 6[l,...,4] 



23 e «- -^S'i i j2 

il...j 4 6[l,...,4] 
hijklmiv) = ViVjVkV{V m 



-2* 



IO / I J1---J5 ' 31 ' 32 'J3 8 J4 8 J5 



:',!2 



ji...i B e[i,...,5] 



- JiL 

;i!2 



+ 



2 3 E 

il...j 5 6[l,...,5] 
hijMmniy) = ViVjVkVlV m V. 



rK rK 
e il.-.J5^ J1 ( 'i J2 i 3 3 lj4lj5 



+ 



ii...j 6 e[i,...,6] 



ii—jeeH.—.e] 



J ~ A K A K A K 

QI03 / j £i'i ...•inOij i- "<.• i Oi. 



ii...j 6 eH,...,6] 



Jl'-'JG i Ji l J2 z J3 l J4 



where we have used the same notation as for the moments. One 
should note that the number of equivalent terms for the moments 
and Hermite polynomials coincides with the factors in the corre- 
sponding terms of the univariate case. The Gram-Charlier coeffi- 
cients can now be calculated by making an ensemble average over 
the Hermite polynomials: 



(40) 

The first coefficient is Co = 1 as in the univariate case and the rest 
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is calculated using the above equation: 





= 


Cij 


= 


( ij k 


i^ijk ' ^ijk 


Cij k I 




Cij k I m 


— f^ijklm 







(41) 



We can find a more compact expression for the last equation 
if we define a skewness term accounting for the asymmetry of the 
distribution function as: 



5M = E K vkhi jh (v) , 



(44) 



ijk 



- y 

3!2 ^ 



Cij klmn 



j 1 ...j 5 £[l,...,5] 
l^ijklm 
f^ijklmn 



« 32 J3 '34 B 



and a kurtosis term accounting for the flatness of the distribution 
function composed by two terms: K, = ICa + ICb with the first 
term given by: 



1 

4!2 



E £j 'l - j e ^31 i 3l *33 *34 Si J 5 »36 

j 1 ...j 6 e[l,...,6] 



+ 



+ 



2 4 3!2 3 J 

^ijklmn 



E 



jl... J6 g[l,...,6] 



E 



31312 C ^l"--'6 K !3 1 l32 l j3 n ' 1 34 l j5'i6 
il...j 6 S[l,...,6] 



3 z 31*32 l J3~ L 34 ? 35 1 36 



4! 

and the second term given by: 
Kniy) = 

1 



KaM = -n^iUjklhijkliv) , 



! j kl 



- y 

6! ^ 



ijklmn 



3!3!2 



E 



J1...J66[1,...,6] 

Now we can rephrase Eq. d43 1 > by: 
P{u) = G(v) [1 + S{y) + K{u) + . . 



(45) 

(46) 

hijklmn (^) ■ 



(47) 



Rearranging terms in the Gram-Charlier series based on the find- 
ings in the univariate case we can build the multivariate Edgeworth 
expansion: 

Au 



ijk 



L E 



720 



P(v) = (det(S))- 1/2 G(iy) 



ijklmn 



(42) 



E 

jl... J6 S[l,...,6] 



ijkl 



Hjklmn 



1 + y^(^»^j'Vfc)c^jfe(^') + -jj ~y j { v i l 'i l 'kVl)chijkl{v) 



1 

6! 



- E 

ijklmn 



1 



3!3!2 



E 



ii...j 8 e[i,...,e] 



Xej'i-ie (^3! ^3-3 M^ j4 ^ J5 ^ 6 )c] if) h i]k i mn {v) + . . . ] 
which can also be written as: 
P(*) = (det(S))- 1/2 G(i^) 

1 + 4 E (*^' $ ^)cE^ /J ^ /2 ^ /2 ^M 



(43) 



i'j'k 1 



ijk 



+ 7T E {*^/*^r)cE^ /2 ^ /2 ^V /2 ^M 



i' j'k'V ijkl 

, E 

i' j' k' V m' n 1 



ji...i 6 e[i,...,6] 



X N S - -1 ^ S ■ ■/ S, , 1 Si,/ ^ S ^ 1 S I h-iikl-mrt 1 ~\~ • • 
/ j it' 23 kk' IV mm' nn' ijretrrtn. \ j 

ijklmn 

where we have inserted the expression for the field u and intro- 
duced the following notation: i' = = l2,k' = i'3,1' = 
«4,m' = i's,n' = i'g. From this expression it is possible to cal- 
culate the probability of a density field <1? given the higher order 
point correlation functions (of the logarithm of the density field!). 



Note that the last term in the previous equation has 10 mem- 
bers (6 indices: 6! /72) which will only be equal under certain sym- 
metry conditions. Our calculations confirm that the first three terms 
in Eq. d47b are a trivial generalization of the univariate case (see 
Eq. |23b . whereas the 4th term is not. 



2.3.1 On the positive definiteness of the expanded Normal 
distribution 

Please note that Pk{v) is not a real distribution function as it is 
not generally positive definite. The truncated expanded distribution 
function should be regarded only as an approximation at fc-th order 
(1st order: lognormal, 2nd order: includes skewness S, 3rd order: 
includes skewness S and kurtosis K,). To ensure a positive definite- 
ness one has to impose the additional condition: 1 + S{v) ^ 0, Vf 
at 2nd order, 1 + S(y) -\-K(y) ^ 0, W at 3rd order, etc.. It is also 
possible to make an exponenzation of the non-Gaussian contribu- 
tion or a quadratic exp ression (see the works on non-Gaussian re- 
alizations in the CMB: IContaldi et ai1l2000l : IContaldi & Magueiiol 
l200ll) : 



P fc (i/) = G{y) exp + K(y) + . 



(48) 



For small skewness and kurtosis terms the Taylor expansion of the 
exponential will be dominated by 1 + S(y) + K,(y) + . . . , so that 
Eq. l |48t will be very close to Eq. \A1\ ensuring positive definite- 
ness. In general one has to define a function T of the non-Gaussian 
contributions which ensures positive definiteness: 



P k [y) = G{v)F{l+S(v)+K{v) + ...) 



(49) 
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Figure 2. The matter statistics depends on the gravitational regime, the cos- 
mic scale and the cosmic time we are looking at. Towards large scales, in 
the linear regime matter is closely Gaussian distributed. At smaller scales 
and later times gravitational clustering will start to depart from Gaussianity. 
When looking at small deviations from Gaussianity a higher order correla- 
tion expansion can be done. The lognormal distribution is a good descrip- 
tion for the further nonlinear regime as it can be regarded as linear in a 
Lagrangian framework which is known to be valid in the quasi-nonlinear 
regime. At even lower scales and late times the lognormal distribution fails 
and an expansion around this distribution function can be done. When shell- 
crossing starts and structures form caustics (the full nonlinear regime) the 
statistical state-of-the-art description is based on numerical N-body simula- 
tions. 



2.3.2 Lognormal limit 

On scales larger than about 10 Mpc the lognormal distribution re- 
sembles verY_well_fte obs£rwdsala}(y statis- 
tics (see lHubblell934l;IWild et all2005l : lKitaura et al.l20"09n . It was 
demonstrated bv lKitaura et alJ ( feOloTT that the lognormal prior can 
be applied at least down to scales of few Mpc to fit the matter statis- 
tics in the overdense regions. Towards the linear regime in the mild 
non-linear regime the correction terms in the Edgeworth expansion 
start to become negligible and we can m odel the density field with 
a multivariate lognormal distribution (see lColes & Joneslll99ll) : 



P(*m|S) 



xcxp 



v /(2 7 r) iv cdet(S) 



(50) 



\ ( ln ( 1 + 5m4 ) ~ ^ si ) S «/ CM 1 + s w) ~ Msj) 



where S is the covariance matrix of the lognormal distribution. Note 



that the covariance matrix is defined by: Sj. 



(SiSj) 



(which does not coincide with the covariance of the overdensity 
field: (SMiS-Mj)) and /x a , describes a constant mean field given by 
(a derivation of both the covariance and the mean can be found in 
appendix let: 



(51) 



with the hats denoting the Fourier transform of the auto-correlation 
matrix. The constant term Su corresponds to the correlation func- 
tion evaluated at zero: S[0], i.e. when it is evaluated for the distance 
of an object at a certain position with itself. 



2.3.3 Gaussian limit 

On very large scales ( > 100 Mpc), in the linear regime, the in- 
formation about the primordial fluctuations has been preserved 
throughout cosmic history. Inflationary scenarios predict a close 
to Gaussian distribution function for the initial density fluctua- 



tions ( seelGuthll98lMGuth & Pill982l:IStarobinskvll982l;lHawkind 
1982 1: iLindd Il982f Ulbrecht & Steinhardtl Il982l : iBardeen et all 
19831) . Therefore we expect to obtain the Gaussian prior in the low 
density limit. When |<5m| <C 1 the signal s can be approximated 
by the linear term of the Mercator series ln(l + 8m j) ~ 5m j and 
the mean is approximately zero ft s , ~ 0. In that case the prior 
distribution for the matter density fi eld is given by the G aussian 
multivariate distribution function (see Bar deen et ail 19861) : 



P(*m|S) 



V(27r)"cdet(S) 



exp 



"2 X/ ^MiSVj^Mj j 
ij / 



(52) 



Please note that the correlation matrix S is now defined by: Sij = 
(^Mi^Mj) ■ This prior is the one required by the Ba yesian approach 
to derive the Wiener-filter (see IZaroubi et alj|l995h a nd has been 
extensively applied t o CMB -mapping (see for example [Burnt et all 
1 1994l : lTegmarld 1997|) and to the large- s cale structure reconst r uction 
(see for example IZaroubi et alj|l995t lErdogdu et al]|2004 120061: 
iKitaura et alj2009l) . 

Intrinsic deviations from the Gaussian or the mild nonlinear 
gravitational regime may be described by expanding the Gaus- 
sian distribution in an Edgeworth expansion. Note that the whole 
multivariate formalism presented in this section can also be ap- 
plied to study weakly no n-Gaussian matter fields as was shown by 
ljuszkiewicz et al.l j 19951) for the univariate case. One has to define 
instead v = S _1//2 <5m and use the higher order correlations cor- 
responding to the overdensity field. However, additional problems 
could arise here yielding negative densities which are not present 
in the lognormal formulation. 

The matter distributions described by the Eqs. d50l > and i52l 
are only valid in a limited range of overdensities. In order to be able 
to go further into the nonlinear regime we need to model the higher 
order statistics. The problem we face here is that we have to intro- 
duce more comple x models increasing the number of parameters 
(see the works b y Scoccimarro et al. 19981 : Taruva & Sodall 19991 : 



(see tne works by acocct 

lMatsubarall2003l : Izhengll2004l: iMatsubar all2008l . including galaxy 
biasing and redshift distortions in higher order correlations). 

The simplest non-trivial higher order correlation representa- 
tion is given by the hierarchical models which we will discuss in 
the next section. 



2.4 Hierarchical model for higher order correlation 
functions 

Hierarchical models rely on the assumption that the cosmological 
structures are to some extent self-similar. In these models higher 
order correlation functions are co nstructed from products of the 
two-point correlation function ( see iFrv & Peebles|[l978l : iFrvll 1984L 
ll986l:lBalian & Schaefferll 19891) : 



6 



•*<n>c=£«»£II& 



(53) 



such that the whole set of points i\ . . , i n is connected by links of 
These links are organized in a tree structure, where a are the 
trees corresponding to each order n. The sum over C a denotes a 
sum over all possible labellings or leafs of a given tree. The re- 
maining freedom is encoded in the hierarchical coefficients for 
each tree a and order n. Note that the hierarchical model always 
refers to the overdensity field 5m, but we are assuming that it also 
applies for the field <1? defined in Eq. (0). 
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Combinatorics of higher order correlation functions in the hierarchical model 



n: order of the correlation function 



m = 2(n — 1): number of indices 



NS 1 : number of trees 



N np = number of leafs (jjjjfc) (i\j\k\l) Nfy (i\j\k\l\m) 

3! /QliHliM _ 4! 2 



(2|1|1)| 3 3=| (3|1|1|1)| 4 4= |}| (4|1[1|1|1)| 5 



(2|2|1|1)| 12 12 = ^2 (3|2|1|1|1)| 60: 



511! 
3! 2 



(2|2|2|1|1)| 60 60 =|^3! 



leafs of the three-point correlation function (2|1[1) | 3 : (£,ij£ik,£ij£jki £ik£jk) 

(2|1|1) (1|2|1) (1|1|2) 

leafs of the first tree of the four-point correlation function (3|1|1|1)| 4 : (£ij£ik£il,Zij£jk£jl,£ik£jk€kl,€iltjl€kl) 

(3|1|1|1) (1|3|1|1) (1|1|3|1) (1|1|1|3) 

leafs of the second tree of the four-point correlation function 

( 2 I 2 |1|1)| 12 : (£y£«k£il>6y&lfjkj£«j£ik£lWj&l£ik^ 

. ' » . ' * » ' * . ' * . ' * . ' 

(2|2|1|1) (2|1|2|1) (2|1|1|2) (1|2|2|1) (1|2|1|2) (1|1|2|2) 

Table 1. The different trees up to order 5 are shown with the corresponding number of leafs/labellings. The trees are constructed by distributing the m indices 
n times so that there is at least one index for each dimension of n neglecting the index order. Note, that by the same procedure one can find that the sixth order 
correlation function has 5 trees instead of 4 as one would naivly expect. To know then the number of leafs we have not only to consider the index order, but 
also the couple combinations of indices. The first factor can be calculated in an analogous way to the number of elements in Fig. {T}. To calculate the latter 
factor one has to first subtract n indices to the tree (since we are interested in the couple combinations we already assume that each correlation function has 
been assigned one index) and then divide (m — n)! by the factorials of the remaining numbers assigned to each index. For example for the second tree of the 
fifth order we would get after subtraction: (2|1|0|0|0), which has three indices left (3! permutations) with 2 equal ones (divided by 2): 3!/2. Dividing this 
number by 3! which comes from the three indices with equal number of appearances, we get that the leafs redundancy number is: JVjL = 2. The leafs for the 
three-point and four-point correlation function are shown in detail at the bottom of the table. 



Particular expressions for the t hree and four- p oint corr elation 
functions were already proposed bv lFrv & Peebles! (see ll978l) . The 
three-point correlation function can then be written as: 



S»x...i 3 — 1 



313233^ [1,2,3] 



(54) 



where we have denoted the only one hierarchical coefficient for the 
three order case as Qz (we show in Tab. Q] how to calculate the 
number of trees and the number of corresponding leafs). Note that 
the three-point correlation function is the sum of the leafs for the 
single tree: 



£sjfc — Qz [£ij£ifc + £ij£jfc + £,ik£jk] 



(55) 



The four-point correlation function has 16 leafs, 4 leafs in the 
first tree and 12 leafs in the second one (see lFrv & Peeblesll 19781 ; 
Bali an & Schaefferll 19891 and our calculation in Tab. [I): 



6i. 



o| ^'i- ••fr&Ji i 32^ i 3X 'is^i' 



Jl-..J4e[l,..-,4] 



+Q4 



il...j 4 S[l,...,4] 



,(56) 



which can also be written as: 

(,ijkl — 

+ Qi [£ij£ik£jl + (.ij^il^jk + ^ij^ik^kl + ^il^ik^kj 
+ ^ik^il^jl + ^ij^jk^kl + ^ik^jk^jl 



(57) 



where we have denoted the first hierarchical coefficient as Q% and 
the second one as Q\. 



2.5 Non-Gaussian multivariate Edgeworth expansion with 
the hierarchical model 

We can use now the three and four-point correlation functions cal- 
culated from the hierarchical model to apply the third order Edge- 
worth expansion to the lognormal probability distribution. Note 
that the results of this section can also be applied for the Gaussian 
case by inserting the correlation functions of the overdensity field 
instead of the correlation functions of the logarithm of the density 
field. In the next subsections we present the skewness and kurtosis 
terms which are required in the Edgeworth expansion (see Eq.|43b. 



2.5.1 Skewness terms 

The second order term in the Edgeworth expansion describes the 
skewness with respect to the lognormal distribution and requires 
the third order Hermite polynomial (see Eq.|39l>. Let us define the 
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weighted Hermite polynomial: 

hi'j'wiy) = J2 S u> /2 S-' /2 S^/ 2 h l]k (v) (58) 



i j k 



with 



7i = E S ^ 



(59) 



Then we can calculate the skewness term 5 as (see Eq. [43] and 
appendix iDl: 

SO) = h^^iikhi jk {v) = h E t,i>j>k'hi>j>k>{v) (60) 



ijk 



i'j'h 1 



3! 



[SiijtSi'k' + Si/ji Sjiy + Si'k , Sj'k']hitj/k'(i') 



i'j'k' 



where we have identified the two-point correlation function 



with Sij . 



2.5.2 Kurtosis terms 

The third order term in the Edgeworth expansion requires the fourth 
and sixth order Hermite polynomials and the three and four-point 
correlation functions. Thus, we separate the kurtosis KZ into two 
contributions KZ = KZa + KZb ■ 

Let us start with the first contribution which requires the fourth 
order Hermite polynomial: 



ijk 

Vt'ilj'Vk'Vi' 



22 E ^-^w a s %f u 

j 1 ...j i e[i,...A] 



- T 

03 



-. . QT-1 tf-1 
03 /, ^]l---3i°i>. i'. J i>. i>. 

A '—f 31 32 33 34. 

31... J4 6[1,...,4] 



The first contribution KZa to the kurtosis term KZ can be written as 
(see Eq.l43l: 

KZa(v) = — ^Kijkihijkiiis) = — ^ £,i>yk>i>hi'j>k>i<{v) ■ 

ijkl i'j'k'V 

which after some calculations (see appendix lElb leads to: 



Qt 
2 



(62) 



+ 



E '/ (i '-- s ' '/, E '/•• S;i .'/. 



(E' 1, ) -'E'" 5 '/ • 2 £ s 



As we saw in sections | |2.2I > and l !2.3t the asymptotic Edge- 
worth expansion has a second term at third order. It is in particular 
the fifth term of the sixth-order moment in Fig. (TJ which is the 
product of two three-point correlation functions. Accordingly, the 
second contribution KZb requires the sixth order weighted Hermite 
polynomial: 

hi'j'k'l'm'n'iv) (63) 

/_ j %v j j kk 11 mm nn' sj/vt/ftftv / 

ijk 



+ 



TT77 E ^ii-ieVi'. Vi'- Vi>. S/i' 

4!2 Z — ' 31 32 33 34 j 5 j 

j 1 ...j 6 e[i,...,6] 
754 E tji—ieVi'. Vi' t / S ir t , 

Z z — . , Jl 32 j 3 j 4 j 5 3q 



ji...j e e[i,...,8] 



QI03 E ^31— J6^i'. i'. ^i'. i'. "S*'. i'. 



ii-..j B e[i,...,8] 



Jl 32 33 34 35 36 



The expression for its inserting the three-point correlation func- 
tion in the hierachical model reads (see Eq.|43t: 



KZb(v) 



- y 



G! 



i' j ' k' I' rn' n' 



3!3!2 



z — ' Jl 32 33 34 35 36 



jl...i 6 6[l,...,6] 



X hi'j'k'l'm'n' i v ) 



— — Q 3 ^ ] [Si'jiSi'k' + Si'j'Sj'k> + Si'k'Sjih'] 

i' j' k 1 V m' n' 

Xhi'j'k>l'm'n'(v) ■ (64) 

After making the corresponding calculations (see appendix IE2b we 
find that the second contribution KZb to the kurtosis term KZ can be 
written as 

K B (S- 1/2 *) = 



(65) 



Qi 



2 

jVi 



-| E ~ 4 E ^iViSjjVi - \ E m®iSijVi$ 

i ij ij 

+J5> S « 1 *3 + §£*? + § (e^ 

+2 .S , //.<!•. + -S , '/,'!' + - ( E 5 '» r ' 

i ij \ i 

+ 7 E ViSijVj — t E ^ — E ^" — s E SuS^S. 



We verify the skewness and kurtosis terms presented here by com- 
paring the univariate case Eqs. ( IA1IIA2I and lA3t to the correspond- 
ing Eqs. dD5l IE 141 and IE37b simplified to a single index (see ap- 
pendices [A][D] and |0. Now the first and second kurtosis contribu- 
tion terms can be added and plugged in Eq. ([43} together with the 
skewness found in the previous section to calculate the probability 
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distribution function of the matter field. We would like to empha- 
size here that the expressions found in this section for the skewness 
and the kurtosis can be efficiently computed (see Egs, 1601 [62] and 
165b . It should be noticed, that they rely on basic operations with the 
fft required to perform convolutions being the most expensive one. 



3 SUMMARY AND CONCLUSIONS 

The avalanche of astronomical data coming from different surveys 
which scan different epochs of the Universe makes it possible to 
map the Universe with unprecedented accuracy. We stress here the 
necessity to model as precise as possible the matter statistics so that 
the multidimensional picture of the Universe can be reconstructed 
with the highest possible fidelity. This implies providing a multi- 
variate higher order statistical description of the structures in the 
Universe. 

We have extended the works done for the univariate case and 
presented the multivariate Edgeworth expansion of the lognormal 
field. We made the expansion explicitly up to third order in pertur- 
bation theory where we had to calculate the multivariate Hermite 
polynomials up to sixth order. The skewness and kurtosis terms in- 
clude the two-point, the three-point, and the four-point correlation 
functions. 

We could show that these terms can be calculated using an- 
alytical expressions for the higher order correlation functions like 
the ones provided by the hierarchical model. 

The expressions derived in this work could be used to gen- 
erate and reconstruct three-dimensional matter fields using higher 
order correlations within a Bayesian framework applying for exam- 
ple the Hybrid Markov Chai n Monte Carlo Hamiltonian sampling 
techn ique (see the works by Jasche & Kitaura 2010; Kitau ra et al.1 
2010). This could have interesting applications to study non- 
Gaussianity in the Large-Scale Structure. 

As new astronomical windows are being opened to map the 
Universe at earlier times in which structures were not clustered as 
much as they are today, we think that the study of the moderate 
nonlinear regime is crucial for an accurate data analysis. Although 
numerical N-body simulations provide a magnificent tool to model 
structure formation, almost without resolution restrictions in com- 
parison to the resolution provided by astronomical observations, 
they do not picture the actual realization of the Universe. There- 
fore, we believe that a special effort should be done to model the 
multivariate statistical nature of the matter distribution which per- 
mits one to extract as much information as possible directly from 
the observational data. We hope that this work serves to contribute 
in this direction. 
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APPENDIX A: UNIVARIATE SKEWNESS AND 
KURTOSIS TERMS IN THE HIERARCHICAL MODEL 



From the hierarchical model ansatz (see section [2T4b we can define 
the univariate three-point correlation function as: £3 = 3Q30- 4 and 
the four-point correlation function by: £4 = 4Qf <r 6 + I2Q4V 6 . Let 
us now look at the skewness and kurtosis terms. 



Al Skewness in the univariate case 

For the skewness we need to define the third order univariate 
weighted Hermite polynomial: /i3(cr _1 $) = a~ z h^{a~ 1 ^) = 
cr~ 6 $ 3 - 3cr~ 4 $. We then get: 

5(<7 ' 1$) = ^ K3hs ( (T ~ 1 ^ = ^3fc3(ff -1 *) = ^(a" 2 c& 3 -3<E-) . 

(Al) 



A2 Kurtosis in the univariate case 

In the case of the kurtosis terms we need to define the 
fourth and sixth order univariate weighted Hermite polynomials: 

hiia- 1 ®) = a^hiia- 1 ^) = cr~ 8 $ 4 - 6cr- 6 $ 2 + 3a' 4 



and h^a' 1 ®) = a^heia- 1 ^) 
45(j- 8 $ 2 - 15cr- 6 . 



15a" iu $ 4 + 



Bl Case k = 1: lognormal 

One finds by shifting the Gaussian integral the solution to the char- 
acteristic function of the lognormal distribution to be given by: 



M*i(t) = (e**)i 



2na 



d$e 2 ^ 



2na 



d$e 2 " 



(B2) 



where we use the following definition: a' 1 = ($ 2 }fe. Recalling the 
expression for <3>: <E> ee log(=) — fi s with p = (p)fcand/z s = { s )k, 
we can write the density as: p = pe* +Ms . Using the characteristic 
function we can calculate all the moments of the density: 



For the lognormal distribution we have: 



(B3) 



<A = 



2%a 



d<P e 2 " 2 p — p e {e )i = p e 



(B4) 



Bl.l Mean 

The particular case for t = 1: 

(p)i = pe" a 

leads to the mean: 



(B5) 



A2. 1 First group of kurtosis terms 

Looking at the fourth order Hermite polynomial contribution we 
get: 



4! 



= \ (Gil + (^ _2 $ 4 - 6$ 2 + 3a 2 ) . 

A2.2 Second group of kurtosis terms 
We then get: 



6! 



6! 



(a- 4 $ 6 - 15a- 2 $ 4 + 45-T 



15a 



(A2) 



(A3) 



B1.2 Variance 
Whereas the case t = 2: 



—2 2^ 3 +2(T 2 



(B6) 



(p)i = p e 

relates the variance of the overdensity 5m to the variance of $: 

„2\ 



(B7) 



- (%)i +1 = e 



p 2 



(B8) 



Thus, one finds that the variance <j 2 of $ under the lognormal as- 
sumption is given by: 



a 2 = ln(<r 2 + 1) , 



wither 2 ee (<&) fc . 



(B9) 



APPENDIX B: UNIVARIATE MEAN AND VARIANCE 

Here we will derive the univariate relations for the mean and vari- 
ance of the expanded lognormal distribution. Note, that the results 
pr esented here are in agreement with the general formula presented 
bv lColombil ( ll994h without an explicit derivation. Let us define the 
fc-th order characteristic function for the variable $ as a function of 
t: 

M* k (t) ee (e'*) fc = J d*fl,(*)e** . (Bl) 

We will look at this expression for the different orders and derive 
from it the corresponding mean and variance. 



B2 Case k — 2: lognormal with skewness 

The characteristic function including skewness yields: 

M* a (t) = (e«*> a 



s/2^a J 



27T<T ./ V 3! 



(BIO) 



1 



2-ko 



d<E> e 



l + i^- 4 (a- 2 $ 3 -3$)') e v 



We have to solve integrals including an exponential term and poly- 
nomials of $. Note however, that we can generate such integrals by 
subsequent derivatives of the characteristic function: 



(e***")! = ^-M*i(t) = -7==- I d$ e -^+**^ 



(Bll) 
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Taking the expression of the generating function for k = 1 (Eq: IC3t 
we can calculate the solutions to the integrals including different 
powers of $: 



B3.1 Mean 



±M n (t) 



toe~ u 

(,1 2 . -,\ 2 tier 2 

(to + 1) o e 2 

(t s o 2 + 3t) o 4 e^' r2 . (B12) 
We can then calculate the characteristic function at order k = 2: 

(B13) 



d 2 , , , v 



M«(t) s =(l + ifs* 8 ! e^*' . 

Defining the integral over the skewness term by: 

S f = ^3t 3 , (B14) 
we can write the second order moment of the density as: 

{ P t )2=fe t ^{e t *) 2 =fe t ^ + t ^ a2 (l + St) . (B15) 

The mean and variance can be obtained by evaluating the latter 
expression at t — 1 and t = 2 respectively. 



B2.1 Mean 



1 2 

^ s = —-o — In (1 + Si) 



(B16) 



B2.2 Variance 

C r 2 =ln[(l + 5 1 )(l+52)" 1 (a| + l)] . (B17) 

B3 Case k = 3: lognormal with skewness and kurtosis 

Here we include also the kurtosis terms: 



M* 3 (t) = <e**) 3 
1 



(B18) 



Since the kurtosis terms use the fourth and sixth order Hermite 
polynomials we need to calculate up to the sixth order derivatives 
of the generating function: 



d* 
dt 4 
d s 
dt 

at 6 



Mi,i(t) = (tV + 6tV 2 + 3)(T 4 e- CT 
-Msi (t) = (tV + lOtV + 15t) o a e^ a2 
M$i (t) = (tV + 15tV + 45tV + 15) cr 6 e^ 



(B19) 



Defining the kurtosis integral by: 



(B20) 



we can write the characteristic function as: 

M*s(t) = (e**) 3 = (1 + S t + Kt) e^" 2 . (B21) 
The t-th order moment of the density yields: 

(p 4 ), = fe^ <e t4, ) 3 = jfe*"'^ " (1+St+Kt). (B22) 

Inserting t = 1 and t = 2 yields the mean and the variance respec- 
tively. 



1 2 

Ms = — -o — ln(l +«Si +/Ci) 



(B23) 



63.2 Variance 

a 2 = In [(1 + 5 a + Ki) (1 + 5 2 + /C 2 ) _1 (<rf + 1)1 . (B24) 



APPENDIX C: MULTIVARIATE MEAN AND 
COVARIANCE 

In this section we derive the multivariate relations for the mean and 
variance of the expanded lognormal distribution. Let us define the 
fc-th order characteristic function for the variable $ as a function 
oft: 

°° 4-91 j.q n 

M* k (ti...t n )= mi---n:) k fr^ (CD 

= (exp^*i$i l j) fc = J d*P fc (*)exp ^E . 

Let us consider now the different fc-order to calculate the corre- 
sponding mean and variance. 



CI Case k = 1: lognormal 

We shift the Gaussian integral to obtain the solution to the charac- 
teristic function of the lognormal distribution: 

M^ih.-.tn) = (exp 

cx Jd* exp f-i *H S hL$im + j 
oc exp f \ tiS ilim t m j /" d* 

V Mm \ I' / \ m' 



x exp 
which simplifies to 



M^{ti . ..t„) = (exp ( j)i = exp ( rE^'i'm'™ 

\ Z / \ Im 

Taking into account that the density field p is related to <1? by: 

Pi = pexp (4>j + p si ), we get: 

(Pn • • • Pi n )i ~P n ex P ( E tl P si i + 2 E tiShimtm J . 



Cl.l Mean 

Setting I = 1 and ti = 1 yields the mean: 

1 



f-l'si 



(C3) 
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CI. 2 Covariance 

Considering now I = 1,2 and ti,ta = 1 we obtain the second 
moment of p: 



1 

\PiPi)l =P ex P ( A^i + Msj + g + ^" + S*i + 



which leads to the covariance: 



Sij = In ({5uiSMj}i + 1) 



(C4) 



(C5) 



C2 Case fc = 2: lognormal with skewness 

Let us recall the expression for the skewness (Eg. 160b: 



(C6) 



(C7) 



ij fc i j fc 

The characteristic function including skewness yields: 
M^ 2 (t!...t n )= <e E < "*'i) a 

7 \ hin, I / 

x 1 + J\ ^Zmhiik(y) ■ 



This integral can be solved by calculating the different moments of 
the lognormal generating function: 



d n 



dt± . . . dt, 



-M^ti.-.tn) oc (C8) 



Performing the derivatives we get: 

^ ^ / ^ 

I \ I'm' 



d 2 



dhdt 2 



MQ^t) = i s ili2 + gi^jt; Si 2 i n 



x exp 



\ V m' 



I I m n 

x exp ( o E tl ' Si i>% n i~ 

\ I'm' 
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This leads to the n-order moment of the density field p: 



(Pn ■ ■■Pi„>2 = p"exp ^ti>n s i v + y^ti'S^, 



+ h E ^ 



(C10) 



x 



which can be simplified to: 



\ i' I'm' / 

X I 1 + — ^ ^ ^iii m i n tltmtn 



C2. 7 Mean 

By setting l,m,n — 1 and ti = 1 we then obtain the mean: 



psi 



C2.2 Covariance 



(CI 2) 



Substituting l,m,n = 1,2 together with ti , £2 = 1 in Eq. dC lib 
we get the second order moment of the density: 



■.'/'</>/) j - /'" <^i> :- s ' , 1(1+ 



(C13) 



which gives us the expression for the covariance: 



= In ((<5Mi<5 M j}2 + 1) + In ( 1 + —^m 



(C14) 



- In ^1 + — (^m + ^jjj + ^uj + + ^jii + ^ijj + ^jij + Cjji, 

Please note, that contracting the expressions found here to the uni- 
variate case gives the same relations as found in the previous sec- 
tion. 



C3 Case fc = 3: lognormal with skewness and kurtosis 

The first contribution to the kurtosis term K, is given by (see 
Eq.HDl: 

£aM = -n^Kijkihijkiiv) = jj ^2&jkihijki(v) , (C15) 



1 J k I 



1 j kl 



6! 
(C9) 



- T 



and the second contribution is 
1 



3!3!2 



3i...j 6 e[i,...,6] 



(CI 6) 



The characteristic function including skewness and kurtosis 
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yields: 

At 43 (ti...t n )=<e E, * , * < «>3 



(C17) 



x (1+S(u)+K(u)) , 

We can then formulate the n-order moment of the density 



field: 

(Pii • 



Pi„}2 = p n exp (y^ty^, + o tl ' Si i' i m' i 



(CI 8) 



+ 6! E 



1 



3!3!2 



E 

.fc 6 6[l,. 



?fcl...fc 6 Ci Jfc 



JkA J/ck -'fcfi 



C3.2 Mean 

The mean is obtained by setting j\ , 



, je = 1 and ti = 1: 



1 



Psi — Sii In ( 1 + ot ^iii + At £,iiii + gj Cii- 



1 

IT 



(C19) 



/ ■ 1 ■ 1 ■ 1 ■ 1 ■ 1 ■ \ 


number of terms 


configurations 




(1 111 lllll) 


l=6!/6! 


^iii^iii 




(1 1 1 1 1 1 1 1 1 2) 


6=6 !/5! 


^iii^iij 1 6— 2 ■ 3 




(1 1 1 1 1 1 2 2) 


15=6!/(4!2) 


£,iii£,ijj 1 6= 2-3i£iij£iij 


9=3-3 


(1 1| 1 Z|Z ZJ 


2U=o!/(J! j ! ) 


siii£,jjj \2^iij^ijj I IS— 


=2-3-3 


(1 1|2|2|2|2) 


15=61/(412) 


6=2-3. 


9=3-3 


(1 2|2|2|2|2) 


6=6 !/5! 


^jjj^jji 1 6=2-3 




(2 2|2|2|2|2) 


1=6 !/6! 







for ji, . 



,36 



1,2 



Table CI. Configurations of j . ; . £4. s 
as needed for the kurtosis contribution to the covariance under certain sym- 
metry conditions (see section IC3.21 . In the first column we consider the 
different configurations of the indices i = i\ and j = 12 disregarding the 
position in the correlation functions. In the second column we calculate the 
number of permutations for each configuration disregarding the position of 
the indices. In the third column we take into account the position of the in- 
dices and identify the different classes of configurations. Note that the sum 
of all the terms gives 64 as we expected. 



The covariance is given by: 

Si] = In ((SuiSyiih + 1) +ln (l + + ^uu + 



(C23) 



C3.2 Covariance 



From Eq. l |C18t we get the second order moment by inserting 
ji,...,je = 1,2 andii,t 2 = 1: 



( r ,,ij .j = p 2 exp (Sij) ( 1 + + lj£ U u + Tir&i ) <C2(» 



1 

IT 



6! 



*{ 1+ 3i Ts+ 4i TK * + ei TK *i 

with the following definitions: 

Ts = £,iii + Cjjj + + + Cjii "t" £sjj + £jij + 

"T" Cjiji Sjjii • (C21) 

Please note that Tk b can be obtained in an analogous way 
from the last term in Eq. dC18t and has 10 x 2 6 = 640 terms. How- 
ever, the number of different classes of terms can be drammatically 
reduced by assuming certain symmetries. In particular assuming 
that any permutation of the indices £»y and the permutation i — s- j 
gives identical terms as it is done in the hierarchical model, leads 
to only 9 classes of terms (see Tab. ICB: f«i|2, £«i£«i|i2=2-2-3, 

(,iii(,ijj 1 12=2-2-3, £,iij 1 18=2-3-3 , £ui£jjj\2, \ 18=2-3-3 With 

all together 64 terms which multiplied by a factor 10 give 640 
terms. The particular expression we find has then the following 
form: 

-^Gb — 20 (2£m -(- 6^iii^jij -\- Q^Hi^ijj ~t~ 9^iij'^iij ~t~ Q^iij^ijj*) • 

(C22) 



APPENDIX D: MULTIVARIATE SKEWNESS TERMS IN 
THE HIERARCHICAL MODEL 

Let us recall the expression for the skewness (Eq. I60t we want to 
calculate here: 

5 (") = y^^ikhijk(v) = ^^2&jkhijk(v) (Dl) 



ijk 



ijk 



y^ [SijSik + SijSjk + Sik.Sjk] hijk(y) 



ijk 



where we have inserted the three-point correlation function from 
the hierarchical model. 

In the following subsection we will calculate the contribution 
of each Hermite polynomial term separately. 



Dl First Hermite term 

Here we go through all the correlation terms applied to the first 
Hermite term: 



Yrijk SijSikViViVk = V, '/.'I'l 

Ei,- fc SijSjkViViVk = E,- m&j 



T,ijk SikSjkVimvk = Efe Vk$l 



S^^iVi, (D2) 



with the factor 3 being due to the fact that the result is identical for 
the three terms of the three-point correlation function. 
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D2 Rest of Hermite terms 

Let us look at the case in which the index of 77 is the same as the 
one which is doubly present in the correlation term SijSik'- 

3x^ SijS ik r)iS]£ = 3x^ S u r)i ■ (D3) 

ijk ii 

This occurs once for each index. Alternatively, the index of 77 coin- 
cides with one of the other two indices: 

3 x J2ij k s ij s jk r ii s jk = 3 x E-y Siji], I „ ^ ^ 

This happens twice. The factor 3 in Eqs. dP3ID4t stands for the 
different leafs/labellings combined with the Hermite terms. 



D3 Result 

Hence, the skewness term is given by: 



5(S" 



-l/2f 



(D5) 



i i i 



APPENDIX E: MULTIVARIATE KURTOSIS TERMS IN 
THE HIERARCHICAL MODEL 

The Edgeworth expansion shows that there are two groups of terms 
contributing to the third order perturbation, which we call the kur- 
tosis: K, = AC a + /Cb- Let us look at each group separately. 



(i) First Hermite term: there are 3 indices for the 3 77's of the 
Hermite polynomial singly coupled to the corresponding index in 
the S functions (in the case chosen here: j, k, I). One index remains 
for the final contraction of the term (here: i): 



4x^ SijSikSariirijTikTii = 4x^ r?i$? 



(E2) 



ijkl 



The factor 4 comes from the fact that this will occur for the 4 in- 
dices which are run by the different leafs of this tree of the four- 
point correlation function. 

(ii) Second Hermite term: 2 indices are assigned to 2 77's singly 
coupled to the corresponding index in the S functions (for example: 
j, k). 2 remaining indices are assigned to 5 1 " 1 with 1 index singly 
coupled to the corresponding index in S (for example: and 1 
index triply coupled (for example: i). This can only happen thrice 
for each term of the four-point correlation function and 4 times for 
each of the indices (4x3): 

4x3x^ SijSikSwntfkSu 1 =4x3x^$?. (E3) 

ijkl i 

The other possibility for this Hermite term is that 1 index of 2 for 
the 2 77's is singly coupled to the corresponding index in S and the 
other one is triply coupled. The other 2 indices are assigned to S~ 
(with occurrence 4x3): 

4 x 3 x y SijSikSu 

ViVjSki = 4x3x^ SuTji^i . (E4) 

ijkl i 

(iii) Third Hermite term: The 2 indices of one of 5* _1 (for ex- 
ample: k, I) are singly coupled to the corresponding indices of the 
S functions. The indices of the remaining S" 1 will be one singly 
coupled (for example: j) and 1 triply coupled (for example: i). This 
happens 3 times for each four-point correlation term: 

4x3x^ SijSikSnS^Ski 1 =4x3xJ]S„. (E5) 

ijkl i 



El First group of kurtosis terms 

The first contribution ICa to the kurtosis term K, is given by (see 
Eq.[43j: 

Ka{v) = 7j Kijkihijkt(f) — 7; ^ i%jkihijki{v) 

ijkl ijkl 

or more specifically: 

/Ca(V) = (El) 

Q4 [SijSikSu + SijSjkSji + SikSjkSki + SuSjiSki] 

+ Q4 [Sij SikSjl + Sij SilSjk ~\~ SijSikSkl ~\~ SilSikSkj 

+SijSuSki + SikSuSji -\- Sij SjkSkl SikSjkSjl 

+SijSjiSki + SaSjkSji + SikSjiSki + SuSjkSki] 
xhijkiif) > 

where we have inserted the four-point correlation function from the 
hierarchical model. 

Since the four-point correlation function has two trees (see 
section [2~4l we will calculate the terms for each tree separately. 

El.l First tree 

The terms corresponding to the first tree are: 
© 0000 RAS, MNRAS 000, 000-000 



El. 2 Second tree 

The terms corresponding to the second tree are: 

(i) First Hermite term: the only possible configuration for each 
four-point correlation term is that 2 77's are singly coupled (for ex- 
ample: k, I) and 2 77's are doubly coupled to the S functions (for 
example: i, j): 

12 x ^ SijSikSjiTurijrikVi = 12 x ^ //,<I>..S'., <!> ;//. , (E6) 

ijkl ij 

thus having a factor 12 for all leafs. 

(ii) Second Hermite term: 2 77's are singly coupled to 2 S func- 
tions (for example with indices: k, I). The remaining indices for 
S^ 1 are doubly coupled (for example: i, j). This happens only once 
for each four-point correlation term: 

12 x Y.^^SjiVkViSrj 1 = 12 x Ij2&i J • (E7) 

ijkl \ i / 

Let us consider now 1 77 is singly coupled (for example: k or /) and 
the other 77 is doubly coupled (for example: i or j) to the S func- 
tions. There are two possibilities: the singly coupled index appears 
in the same S function with the doubly coupled index: 

12 x 2 x S l] S rk S J ir ll r,kS- l 1 = 12 x 2 x ^ Si^i , (E8) 

ijkl ij 



18 Kitaura, Francisco-Shu 



or the singly coupled index appears in the remaining S function: 



E2 Second group of kurtosis terms 



12 x 2 x J2 S lJ S lk S J ir lt ViS~ L 1 = 12 x 2 x Y Su^i . (E9) 

ijkl ij 

There are 2 configurations for both cases for each four-point corre- 
lation term, hence 12 x 2. 

The last configuration for the second Hermite term is based on 2 
indices of r\ doubly coupled to the S functions (this happens only 
once for each four-point correlation term: 12 times in total): 

12 x ^ SijSiuSjmriiSu 1 = 12 x Y , (eio) 

ijkl ij 

(iii) Third Hermite term: the last Hermite term has two possibil- 
ities: either only one of the indices of the S^ 1 functions coincides 
respectively with one of the indices of the S functions (this can 
only happen once for each four-point correlation term): 

12 x Y S t3 S tk SjiS- l 1 S- k 1 = 12 x Y s " , ( E1 !) 

ijkl i 

or two indices of one of the 5 1-1 functions coincide with the corre- 
sponding indices of one of the S functions (this can occur twice for 
each four-point correlation term: 12 x 2): 

12 x 2 x Y SijSikSjiS^Su 1 = 12 x 2 x Y s v > ( E12 > 

ij k I ij 



El. 3 Result 

Putting the contribution of both trees together we get the first kur- 
tosis term: 



£ A (S- 1/2 *) = 



(El 3) 



91 

2 



+ os) £ s »$wi + \ (qi + oi) £ s - 



+ 



£ Vi®iSij®jVj - £ '/ s 



2 

ijVj 



£ *• -V'" - s '/ • 2 £ s ^' 

i / ij ij 

or equivalently: 

£a(S- 1/2 *) = 

^E^E^-E^ 



(E14) 



2 



^ + Q") £ E 5 ^ + \ («4 + Q 4 b ) E s - 



+■ 



££s« 1 # i *<S y * i £sr/* j 



EE^ : «K.v'£.v >. 

ij i j 

(£"••) - 2 E •" • s E V*i + 2 E s » 



The second contribution /Cb to the kurtosis term K, is given by (see 
Eq.gUi: 



1 



- y 



(i! 



(El 5) 



3!3!2 



n...j a 6[i,...,6] 

xAil...i 6 (f) 

— "^|"Q3 ^ ^ ['S'ijSifc + SijSjk + 5'jfc5'jfc] 
ijfcimn 

X ['S'zmS'zn + Sl m S mn -(- »S/ n 5' mn ] 
X foijklmn (^) 5 

where we have inserted the three-point correlation function from 
the hierarchical model. 

The three-point correlation function which has only one tree 
in the hierarchical model and the sixth order Hermite polynomial 
which has four terms. Let us partition the problem into these four 
Hermite terms taking into account the symmetries intrinsic to the 
hierarchical three-point correlation function. 

E2. 1 First Hermite term 

Here we always have 4 77's singly coupled and 2 77's doubly coupled 
to S functions. Since we have the three-point correlation function 
squared we will have 3x3 terms: 

9 x 

ijklmn 

2 

(E16) 



= 9x (X>**) 



E2.2 Second Hermite term 

(i) 4 jy's are singly coupled to the S functions leave the 2 remain- 
ing indices for the function doubly coupled to the S functions 
(this can happen only once for each of the 9 hierarchical correlation 
terms): 

9 X £ SijSikSlmSlnTljrikrjmrinS^ 1 
ijklmn 



f> E <i, "- S ', l<i> ? 



(El 7) 



(ii) 3 77's are singly coupled to the S functions. The remaining 
77 is doubly coupled to the S functions. It can happen that the re- 
maining 77 is coupled to the S functions which are also coupled to 
2 77's: 

9 x 4 x £ SijSihSi m Si n rjjr]kr] m rnS^ 

ijklmn 



i j 

or only to 1 : 

9 x 4 x Y, SijSikSimSiniljr)kr)mriiS~^ 

ijklmn 

= 9x4x^ $f rji . 



(El 8) 



(El 9) 
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(iii) 2 77's are singly coupled to the S functions. Both singly cou- 
pled 77 indices are coupled to S functions which have a common 
index (can happen twice): 



9 x 2 x ^2 S i:j SikSi m Si n riir] j ri k r]iS rrl 1 n 



ijklmn 



Both singly coupled 77 indices are coupled to S functions with dif- 
ferent indices (has 4 combinations): 

9 X 4 X ^2 SijSikSl m Sl„riiT)jTl m T)lSkn 
ijklmn 

= 9 x 4 x 'I 'I'.^ , '/, '!', • (E20) 



E2.3 Third Hermite term 

(i) 2 tj's are singly coupled to the S functions. Both singly cou- 
pled 77 indices are coupled to S functions which have a common 
index: 

9 X 2 X ^ y Sij SikSlmSln^lj^kSn S mn 
ijklmn 

= 9x2x^S lt ^S^, (E21) 

i j 

or by a permutation of the indices of the S^ 1 -functions: 

9 X 4 X ^ ^ Sij SikSlmSlnTljT)kS^ rn Si rl 
ijklmn 

F.2 



9x4x^]$? 



(E22) 



(ii) 2 77's are singly coupled to the S functions. Both singly cou- 
pled 77 indices are coupled to S functions which have no common 
index: 

9 X 4 X ^ ^ SijSikSlmSlnVjVmS^i^, Si n 



ijklmn 



9x4x (l>^ > 



(E23) 



or by a permutation of the indices of the S 1 -functions: 

9 X 8 X ^ ] SijSikSlmSlnVjVmS il S kn 



ijklmn 



= 9 x 8 x $ ^ 



(E24) 



(iii) 1 77 is singly coupled to 1 S function. 1 77 is doubly coupled 
to 2 S functions with one of them sharing the same index as the 
singly coupled one. The rest of the indices which appear only once 
and do not share the same doubly present index in the S functions 
do not mix in the S^ 1 functions (4 combinations in the 77 indices): 



9x4x ^2 S i jSikSi rn Si n riir]jS kl 1 S m l l 

ijklmn 
= 9 X 4 X SiiVi^i ■ 



(E25) 



The rest of the indices which appear only once and do not share the 
same doubly present index in the S functions are mixed in the S^ 1 



functions (4 combinations in the 77 indices and 2 combinations in 
the S^ 1 indices): 

9 x 8 x SijSikSimSinriirijS km S ln 

ijklmn 



= 9 x 8 x 



(E26) 



(iv) 1 77 is singly coupled to 1 S function. 1 77 is doubly coupled 
to 2 S functions with non of them sharing the same index as the 
singly coupled one. The rest of the indices which appear only once 
and do not share the same doubly present index in the S functions 
do not mix in the S^ 1 functions (4 combinations in the 77 indices): 



9 x 4 x ^2 SijS ik SimSiniljriiS ik S mn 

ijklmn 

= 9 x 4 x ^2 Surji^i . 



(E27) 



The rest of the indices which appear only once and do not share 
the same doubly present index in the S functions are mixed in the 
S^ 1 functions (4 combinations in the 77 indices and 2 combinations 
in the S^ 1 indices): 

9 X 8 X ^2 SijSikSlmSlnVjrilSrmSkn 
ijklmn 



9 x 8 x ^2 S ll r) i $ l . 



(E28) 



(v) Both 77's are doubly coupled to the S functions. The rest of 
the indices which appear only once and do not share the same dou- 
bly present index in the S functions do not mix in the S~ 1 functions 
(only 1 combination): 

9 x ^2 SijSikSimSi n r]iriiSj k S mn 

ijklmn 



= 9 x {^2Sur,}j 



(E29) 



The rest of the indices which appear only once and do not share the 
same doubly present index in the S functions are mixed in the S^ 1 
functions (2 combinations in the S^ 1 indices): 

9 X 2 X ^ ^ SijSi k SlmSln r QiT)lSj rn S kn 
ijklmn 

= 9x2x53^5?-%. (E30) 



E2.4 Fourth Hermite term 

(i) Both indices of 2 functions are pairwise the same as the 
indices of 2 5 functions (4 combinations): 

9 X 4 X ^ ^ Sij SikSlmSlnS^j Si m S kn 
ijklmn 



(E31) 



(ii) Both indices of 1 S 1 function are the same as the indices 
of 1 S function (4 combinations): 



9 X 4 X ^ ^ Sij SikSlmSlnS^j S k i S mn 
ijklmn 

= 9 x 4 x J2 S " ■ 



(E32) 
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(iii) Non of the indices of the S^ 1 functions coincides pairwise 
with the indices of the S functions. The indices which are doubly 
present in the S functions do not coincide with the indices of a 
single S" 1 function: 



9x4 E] Sij SikSi m Si n S ir ^ S S^ 1 

ijklmn 



(E33) 



(iv) Non of the indices of the S^ 1 functions coincides pairwise 
with the indices of the S functions. Both indices which are doubly 
present in the S functions coincide with the indices of a single S~ 
function. The rest of the indices which appear only once and do not 
share the same doubly present index in the S functions do not mix 
in the S" 1 functions (only 1 possibility): 



9 X S ^ Sij SikSlmSlnS^i Sjfc Smn 
ijklmn 

= 9 X y ^ SiiS i; j Sjj . 



(E34) 



The rest of the indices which appear only once and do not share 
the same doubly present index in the S functions are mixed in the 
S^ 1 functions (2 combinations): 



9 x 2 x SijSikSi m Si n S a Sj m S kn 

ijklmn 



which can also be written as: 

/Cb(S" 1/2 *) = 



(E37) 



Ql 



V i j / ij 

ij j i j 

-iE^E^^E^ 

ij i j 

-\ E E srM&iSij E •*>'.>."•• 

ij i j 

+iE^^ 2 + ^E^ + ^(E^) 2 

ij i \ i / 

+2 e s u e .s + e e s»*j*i 

i j ij j 

+1 (E su E 2 + i E E E 

\ i j / ij i j 

— J E ~ E S*> ~ o E SiiS *j 1S 3- 



Please note that one gets the same result as in Eqs. JA1IIA2I 
and lA3l l by simplifying the corresponding Eqs. JD5IIE141 and lE37b 
to a single index. This actually demonstrates that we have gone 
through all the possible configurations of indices as the number of 
equivalent Hermite terms is recovered. 



(E35) 



E2.5 Result 

Summing up the terms we get: 

£b(S" 1/2 *) = (E36) 



Ql 



5 I)*?* -gE*W«5-5E*«? 



~2 E <1 ''''' - 4 E <1 '-'/ - ..'/. - \ E Vi$iSijVi®i 

i ij ij 

+2^&'(,*, + E s + g ( E ) 

i ij \ i / 

- E ViSijVj — j E ^ — E ^" — » E ^"^j ^j. 



4^" 1JU 4 
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